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Abstract  1/^/ 

We  present  a  new  approach  towards  the  construction  of  a  genuinely  multidimensional 
high-resolution  scheme  for  computing  steady-state  solutions  of  the  Euler  equations  of  gas 
dynamics.  The  unique  advantage  of  this  approach  is  that  the  Gauss-Seidel  relaxation  is  stable 
when  applied  directl}^  to  the  high-resolution  discrete  equations,  thus  allowing  us  to  construct 
a  very  efficient  and  simple  multigrid  steady-state  solver.  This  is  the  only  high-resolution 
scheme  known  to  us  that  has  this  property.  The  two-dimensional  scheme  is  presented  in 
detail.  It  is  formulated  on  triangular  (structured  and  unstructured)  meshes  and  can  be 
interpreted  as  a  genuinely  two-dimensional  extension  of  the  Roe  scheme.  The  quality  of  the 
solutions  obtained  using  this  scheme  and  the  performance  of  the  multigrid  algorithm  are 
illustrated  by  the  numerical  experiments.  Construction  of  the  three-dimensional  scheme  is 
outlined  briefly  as  well. 
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was  also  supported  in  part  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract 
No.  NASl-19480  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications  in  Science 
and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681. 
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1.  Introduction*  A  need  for  a  fast  and  accurate  steady-state  compressible  flow 
solver  exists  in  many  areas  of  science  and  engineering.  Such  a  need  is  partictdarly 
acute  for  the  problem  of  aerodynamic  design.  In  this  case  steady-state  solutions  of 
the  compressible  flow  past  bodies  have  to  be  computed  repeatedly,  each  time  with 
variations  in  the  body’s  geometry.  This  allows  us  to  find  the  shape  with  optimal 
aerodynamic  parameters  relying  on  computations  only.  Thus  the  necessity  of  costly 
wind-tunnel  experiments  can  be  largely  reduced. 

The  search  for  genuinely  midtidimensional  schemes  for  the  compressible  Euler 
equations  was  motivated  by  the  expectation  that  they  will  provide  new  possibilities 
(compared  to  the  dimensionally  split  schemes  in  common  use  now),  like: 

•  capturing  physics  of  the  fluid  flow  more  accurately; 

•  constructing  a  more  efficient  steady-state  (multigrid)  solver. 

Considerations  regarding  the  first  point  can  be  found  in  [14], [15].  The  main  motivation 
for  constructing  the  truly  multidimensional  scheme  in  this  work  is  the  improvement 
of  multigrid  efficiency.  It  was  observed  in  [20]  that  pointwise  Gauss-Seidel  relaxation 
is  unstable  when  applied  directly  to  the  high-resolution  dimensionally  split  scheme 
even  in  the  simple  case  of  a  two-dimensional  scalar  advection  equation.  Therefore,  the 
steady-state  Euler  solver  constructed  in  [21]  relies  on  the  defect  correction  technique, 
which  is  not  a  fully  efficient  way  to  utilize  multigrid  methods.  Another  possibility 
to  avoid  this  difficulty  is  to  use  a  multigrid  algorithm  that  employs  a  well-known 
midti-stage  Runge-Kutta  relaxation  technique  which  was  developed  in  [9], [7], [8]. 

However,  any  further  improvement  of  the  multigrid  efficiency  requires  us  to  address 
this  problem  directly:  it  is  necessary  to  develop  a  new  high-resolution  (at  least  at 
the  steady-state)  discrete  scheme,  such  that  Gauss-Seidel  relaxation  is  stable  when 
applied  directly  to  the  resulting  discrete  equations.  This  wais  the  main  motivation  for 
constructing  the  genuinely  two-dimensional  advection  scheme  in  the  control  volume 
context  in  [17], [18].  The  novel  feature  of  this  scheme  was  the  two-dimensional  limiter, 
i.e.  the  limiter  function  that  relies  on  the  ratio  of  two  finite  differences  in  different 
directions. 

There  has  also  been  developed  another  class  of  genuinely  two-dimensional  ad¬ 
vection  schemes  -  the  so-caUed  ‘ffiuctuation-splitting”  schemes,  (see  [6],[22]).  These 
schemes  are  also  equipped  with  a  certain  nonlinear  mechanism  that  allows  to  combine 
high-resolution  and  a  positivity  properties.  Several  variants  of  such  a  nonlinear  mech¬ 
anism  were  devised  using  some  geometric  considerations.  The  remarkable  feature  of 
this  approach  is  the  simplicity  of  the  schemes  formulated  on  unstructured  triangidar 
grids. 

The  strong  relationship  between  the  fluctuation-splitting  and  control  volume  ad¬ 
vection  schemes  was  established  in  [19].  As  a  result,  the  fluctuation-splitting  scheme 
utilizing  two-dimensional  limiters  was  constructed.  The  action  of  such  a  limiter  func¬ 
tion  can  be  given  a  purely  algebraic  interpretation. 

Even  though  some  genuinely  two-dimensional  advection  schemes  were  available 
already  several  years  ago,  the  task  of  extending  these  ideas  to  the  systems  of  equations 
appeared  to  be  rather  difficult. 

The  ‘‘algebraization”  of  the  advection  scheme  formulation  appeared  to  be  crucial 
for  the  purpose  of  this  work  -  constructing  a  truly  multidimensional  scheme  for  the 
Euler  system.  The  resulting  scheme  is  capable  of  producing  high-resolution  steady- 
state  solutions.  Its  unique  feature  is  that  the  Gauss-Seidel  relaxation  is  stable  when 
applied  directly  to  the  high  resolution  discrete  equations. 
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The  paper  is  organized  as  follows:  in  §2  we  give  a  brief  introduction  into  the 
fluctuation-splitting  approach  for  the  scalar  advection  equations  and  present  consid¬ 
erations  for  extending  this  approach  to  systems  of  equations  in  two  dimensions.  In 
§3  we  take  a  closer  look  at  the  Roe  scheme  for  the  Euler  equations  in  one  dimension. 
Then  we  construct  a  truly  two-dimensional  scheme  for  the  Euler  systems  on  structured 
(in  §4)  and  unstructured  (§5)  triangular  grids.  In  §6  we  outline  first  as  a  preliminary 
the  construction  of  the  three-dimensional  advection  scheme.  Then  we  present  a  truly 
three-dimensional  scheme  for  the  Euler  system.  §7  describes  the  multigrid  algorithm 
that  employs  the  constructed  truly  two-dimensional  scheme.  Numerical  experiments 
are  presented  in  §8,  followed  by  discussion  and  conclusions  in  §9. 

2.  Preliminaries  and  motivation.  In  the  first  part  of  this  section  we  illus¬ 
trate  the  fluctuation-splitting  approach,  first  on  the  example  of  scalar  advection  in 
one  dimension.  Then  we  present  the  construction  of  a  simple  truly  two-dimensional 
advection  scheme  that  relies  on  two-dimensional  limiters.  In  the  second  part  we  dis¬ 
cuss  the  difficulties  that  arise  when  trying  to  apply  the  scalar  advection  schemes  to 
discretize  systems  of  equations  in  multidimensions. 

2.1.  Fluctuation-splitting  approach.  We  shall  give  here  a  brief  description 
of  the  fluctuation-splitting  advection  scheme  in  one  and  two  dimensions  (for  details 
see  [6], [19]).  These  schemes  are  required  to  have  properties  of  positivity  and  linearity 
preserving.  We  shall  define  here  the  positivity  property. 

Definition  2.1.  A  scheme  is  said  to  be  of  the  positive  type  if  any  solution  value 
on  the  new  time  level  obtained  by  this  scheme  can  be  written  as  a  positive  combination 
of  the  values  from  the  previous  time  level. 

Solutions  obtained  by  using  positive  schemes  wiD  satisfy  a  certain  maximum  prin¬ 
ciple  and,  therefore,  do  not  exhibit  oscillatory  behavior  in  presence  of  discontinuities. 

The  notion  of  linearity  preserving  is  used  to  characterize  high-resolution  at  the 
steady-state  schemes.  It  is  trivial  in  one  dimension.  Therefore,  we  shall  give  its  precise 
definition  in  §2.1.2. 

2.1.1.  Advection  in  one  dimension.  Consider  a  linear  one-dimensional  ad¬ 
vection  equation 

(1)  ut  +  auj:  =  0 

We  are  interested  in  solving  it  numerically  on  a  grid  with  meshsize  h  (see  Fig.l). 
Fluctuation  is  defined  as  the  residual  of  the  equation  on  the  linear  element  multiplied 
by  the  volume  of  the  element.  In  our  case  the  fluctuation  on  the  segment  [i  —  l,i]  is 
defined  as  the  residual  of  (1)  on  this  segment  multiplied  by  its  length  h 

R  =  --a{ui  —  Ui-i) 

The  numerical  solution  can  be  interpreted  as  a  two-stage  process 

•  compute  the  fluctuation  on  each  element  and  distribute  (split)  it  among  the 
vertices  of  this  element; 

•  update  the  solution  at  each  vertex  due  to  the  accumulated  portions  of  fluctu¬ 
ations. 

Denoting  by  r  the  time-step  and  distributing  the  fluctuation  on  the  segment  [i  —  l,i] 
equally  between  the  adjacent  gridpoints 

hu^+^  +'^R  +  C.O.E. 
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hu^+'  =  huf  +^R  +  C.O.E. 

we  obtain  the  central  scheme,  which  is  neither  positive  nor  stable,  though  second  order 
accurate  in  space.  Here  “C.O.E”  stands  for  “contributions  from  other  elements”  and 
will  be  omitted  in  the  remainder  of  this  paper. 

The  upwind  scheme  can  be  obtained  by  adding  a  proper  amount  of  the  artificial 
viscosity  to  the  central  scheme.  The  fluctuation  distribution  formulae  thus  become 

(2)  hu:+^  =  huU  +^{R-R) 

where 

5  =  — |c|(Mi  -  =  sign(a)E 


is  the  artificial  viscosity. 

It  is  easy  to  see  that  in  this  case  the  entire  fluctuation  on  the  segment  [i  —  l,i] 
contributes  to  the  solution  update  either  at  the  left  or  at  the  right  nodes  of  the  segment 
depending  on  the  advection  direction  (sign). 

2.1.2.  Two-dimensional  advection.  Consider  a  linear  two-dimensional  advec¬ 
tion  equation 


Ut  -f  aux  +  buy  =  0 

The  fluctuation  on  the  triangle  T  (see  Fig.2)  is  given  by 

(3)  R=R^  +  Ry, 

where 


R""  =  -|K«o  -  «3)] 

Ry  =  -|[6(w3  -  U4)] 

Assuming  that  the  fluctuation  distribution  formulae  in  this  case  are 

+  ^(R^  +  R^) 

(4)  =  h‘^u^  +  I [(i?®  -  R^)  +  (Ry  +  Rv)] 

^  f^2yn  ^ 

where  R^,Ry  are  the  artificial  viscosity  terms  defined  by 

/c\  R^  =  siga(a)R^ 

Ry  =  sigTL{b)Ry, 

we  obtain  the  dimensional  upwind  scheme  which  is  positive,  but  only  first  order  accu¬ 
rate. 

Definition  2.2.  The  fluctuation-splitting  scheme  is  called  linearity  preserving  if 
whenever  the  fluctuation  on  the  triangle  T  vanishes  then  the  scheme  leads  to  a  zero 
update  in  each  of  the  three  vertices  of  the  triangle. 

It  was  observed  in  [6]  that  linearity-preserving  schemes  on  structured  grids  produce 
numerical  solutions  that  are  second-order  accurate  in  the  steady  state. 
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Introduce  the  following  quantities 


(6) 


=  R=^  +  Ry^iQ) 
Ry'  =  Ry  +  R^^ 


where 


(7) 


Q  =  - 


R^ 

Ry 


and  $  is  a  Lipschitz  continuous  limiter  function  such  that 


(8)  0  <  ^(Q)  <1,  0  <  <  1 

and 


(9)  'f(l)=l 

Substituting  for  R^^Ry  into  (4)  and  (5)  we  obtain  a  linearity  preserving 

(second  order  accurate  at  the  steady-state  in  the  case  of  structured  grid)  scheme. 
Using  the  following  identity 

(10)  = 

it  is  easy  to  see  that  the  scheme  defined  by  (4), (5)  and  (6)  is  of  the  positive  type. 

It  is  also  obvious  from  (10)  that  such  scheme  is  conservative  because 

R^^  +Ry*  =R^  +  Ry  =  R 


(for  more  details  see  [19]). 

2.2.  Hyperbolic  systems  of  equations.  We  shall  explain  here  what  is  the 
main  obstacle  encountered  when  constructing  a  truly  two-dimensional  numerical  scheme 
for  a  hyperbolic  system  of  equations.  We  shaU  also  describe  our  approach  to  overcom¬ 
ing  it.  An  understanding  of  the  basic  differences  between  the  one-  and  two-dimensional 
cases  is  required  for  this  purpose.  These  differences  will  be  illustrated  on  linear  systems 
of  equations. 

2.2.1.  One  dimension.  Consider  a  system  of  conservation  laws  in  one  dimen¬ 
sion 

(11)  ut  +  f{u)j,  =  0 

For  the  purpose  of  the  discussion  here  it  is  sufficient  to  look  at  the  linear  constant 
coefficient  case  of  (11).  Consider  the  following  system  of  partial  differential  equations 

(12)  ut  +  Auj,  =  0 

where  A  is  an  A  x  A  matrix.  The  system  (12)  is  said  to  be  hyperbolic  if  matrix  A 
has  a  complete  set  of  real  eigenvalues. 

Denote  by  T  the  matrix  of  right  eigenvectors  of  A,  Then 

A  =  T-MT 
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is  a  diagonal  matrix.  Introducing  characteristic  variables 

w  =  T~^u 

(12)  can  be  rewritten  as  a  set  of  N  decoupled  advection  equations: 

(13)  Wt  +  Awa:  =  0. 

It  is  clear  from  (13)  that  a  one-dimensional  advection  scheme  can  be  applied  in  a 
straightforward  way  to  solve  a  (linear)  hyperbolic  system  of  equations  in  one  dimen¬ 
sion.  A  discussion  concerning  the  nonlinear  Euler  system  in  one  dimension  will  be 
presented  in  the  §3. 

2. 2. 2*  Two  dimensions.  A  linear  system  of  partial  differential  equations  in  two 
dimensions  of  the  following  form 

(14)  Ut  +  Aua:  -h  Buy  =  0 
is  said  to  be  hyperbolic  if  the  matrix 

A  =  cos  (t>A  +  sin  4)B 

has  a  complete  set  of  real  eigenvalues  for  :  0  <  <;^>  <  27r.  In  this  case  there  exist 
matrices  Ta  Q-nd  Tq  such  that  AT  a  2ind  T^^BT^  are  diagonal.  This  is  usually 

utilized  to  construct  the  so-called  dimensionally  split  schemes. 

However,  matrices  A  and  B  in  general  do  not  commute.  Therefore,  Ta  Tb-)  i.e. 
a  hyperbolic  system  in  two  dimensions  cannot  be  represented  as  a  set  of  decoupled 
advection  equations  (unlike  the  one-dimensional  case).  This  means  that  a  truly  two- 
dimensional  advection  scheme  cannot  be  applied  in  a  straightforward  way  to  discretize 
a  hyperbolic  systems  of  equations. 

Much  of  the  research  effort  in  the  last  several  years  concentrated  on  finding  a  way 
to  apply  the  multidimensional  advection  schemes  to  discretize  hyperbolic  systems  of 
equations,  in  particular  the  Euler  equations  of  gas  dynamics.  One  of  the  major  direc¬ 
tions  was  the  so-called  wave  modeling  (see  [14], [5]).  This  approach  concentrated  on 
finding  a  way  to  represent  (locally)  the  physics  of  two-dimensional  flow  of  a  compress¬ 
ible  fluid  by  a  finite  number  of  simple  waves,  each  one  having  an  associated  advection 
equation. 

The  approach  we  present  here  is  concerned  not  with  applying  an  advection  scheme 
to  discretize  a  system  of  equations  in  two  dimensions,  but  rather  with  applying  to  the 
systems  of  equations  the  same  strategy  that  was  used  when  constructing  a  scalar 
advection  scheme.  The  construction  of  the  genuinely  two-dimensional  scheme  for 
the  Euler  will  be  presented  in  §4  for  the  structured  grids  and  in  §5  for  the  general 
unstructured  grids  case, 

3.  Euler  equations  in  one  dimension  and  Roe  scheme.  The  Euler  equa¬ 
tions  of  gas  dynamics  in  one  dimension  can  be  written  as  follows 

(15)  Ut  +  F{u)^  =  o, 
where 

(  pu  \ 

;  F{u)  =  pv?-^p  \, 

\  puH  ) 
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the  enthalpy  H  is  defined  by 


(1.7) 

the  speed  of  sound 


H  = 


e  +  p 


- 1-  -■ 

7-12  ’ 


and  the  pressure 

(19)  P  =  (7- l)(e- ^pw^) 


The  Jacobian  matrix 


(20) 


A  = 


dF 

du 


has  real  eigenvalues  only. 

In  order  to  apply  the  upwind  differencing  approach  to  the  nonlinear  system  (15) 
a  particular  conservative  linearization  procedure  was  introduced  by  Roe  [13].  For  two 
states  ui  and  a  matrix  A{ui^Ur)  was  constructed  having  certain  properties  called 
together  Property  U,  We  address  the  reader  to  [13]  for  the  details.  Let  us  only  mention 
here  one  of  these  properties:  the  following  identity  holds  for  any  ui^Ur' 


(21) 


A{uuUr)  •  {Ur  -  Ul)  =  F{Ur)  - 


3.1.  Fluctuation-splitting  formulation.  Define  the  fluctuation  on  the  consid¬ 
ered  element  (segment  [t  —  l,t],  see  Fig.l) 


(22) 


Following  the  Roe-linearization  procedure  we  assume  that  the  quantities  which  vary 
linearly  on  the  segment  [z  —  l,i]  are  the  elements  of  the  so-called  parameter  vector 
(see  [13]) 

(23)  m  =  {mi,m2,Tn3)'^  =  y/p{l,u,H)'^. 


Therefore,  its  average  value  on  [z  —  1 ,  z]  is 

mj_i  -|-  m, 

2 

Thus  the  following  Roe-averaged  quantities  can  be  defined 


p*  =  ml 

(25)  u  =  fh^lrhi 

H  =  ms/rhi 


Introducing  the  Roe-averaged  Jacobian  A  =  A(ui-i,Ui)  and  recalling  (21)  we  can 
rewrite  the  fluctuation  (22)  in  the  following  form 


(27)  R  =  -A  •  (ui  -  Ui-i). 

Having  in  mind  the  construction  of  an  upwind  scheme,  it  is  convenient  to  use  the 
following  representation  of  the  matrix  A  (see  [13]) 


where  A  is  a  diagonal  matrix 


where 


A  =  TAf-^, 


/  Ai  0  0 

A  =  0  Aa  0 

I  0  0  A3 


Ai  =  U,  A2,3  =  u  ±  c 


1  2  3  ^ 

are  the  eigenvalues  and  T  =  {E  ,E  ,E  )  is  the  matrix  of  right  eigenvectors  of  A 


E^=(  u  |;  E^=(  u-c  ■  E^^i  u  +  c 

\  S^/2  /  \  H  —  uc  I  \  H  +  uc 


The  fluctuation-splitting  formulation  of  the  Roe  scheme  can  be  written  as  follows 

^31)  h<_+i  =  h<_i  +^(R+R) 

^  =  hvl>  -bf  (iJ  -  R), 


where 


is  the  artificial  viscosity,  and 


R  =  -f  |A|(tni  -  t»,_i) 


w  =  T 


are  the  so-called  characteristic  variables. 

Alternatively,  the  artificial  viscosity  can  be  also  written  as 

(34)  R  =  -\A\{ui-Ui.i) 


R  =  — sign(^)J2 
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3.2.  Roe  scheme  revisited.  It  is  convenient  for  the  purpose  of  this  work  to 
introduce  the  foUowing  auxiliary  variables  where 

(36)  ds  =  dp - ^ 

is  the  entropy.  The  non- conservative  formulation  of  the  Euler  equations  in  these 
variables  takes  the  following  form 


(37) 


St  +  US^  =  0 

pUt  +  pUUj:  +  =  0 

Pt  +  UPa:  +  pC^Uj:  =  0 


Introduce  the  fluctuations  of  the  Euler  equations  in  the  auxiliary  variables 
(38) 


Si  -  Si-i 

r  =  -A\ 

P,  -  Pi-1 


where 


(39) 


It  can  be  easily  verified  that 

(40) 

where 


R  =  CaT, 


(41) 


Ca  = 


^  1  0  l/c^ 

u  1  u/c^ 

^  w^/2  u  1/(7  -  1)  +  ti^/(2c^) 


Therefore,  the  Roe  scheme  (31)  can  be  written  as  follows 
(42) 


=  hu^_-i  +^Ca(r  +  r) 
=  hu'l  +|Ca(r-r) 


where  the  artificial  viscosity 


(43) 
with 

(44)  = 
We  can  also  write 

(45) 


f  =  C-^R  =  -C-'r|A|K  - 


l-(7-l)t2/(2e^)  (7-l)«/c^  -(t  -  1)/^^ 
-u  1  0 

(7- I)it2/(2c2)  _(7-l)u  7-1 


Si  ^1  —  1 

/>*(««•- “i-i) 
Pi  -  Pi-1 
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or 


r  =  sigii(j4)r 


It  is  evident  from  (43)  and  (45)  that  the  matrix  of  right  eigenvectors  of  A  is 


err 


or  T  =  (ei,e2,e3)  where 


ei  =  0 

V  0 


€2  =  -C 


63=  C 


Remark  3.1.  It  is  obvious  from  the  structure  of  the  matrix  T  that  the  auxiliary 
variable  s  coincides  with  the  first  element  of  the  characteristic  variable  vector  (33) 


Denote 


wi  =  s 


M  =  sign(J[) 


In  order  to  make  a  computer  program  more  efficient  it  is  useful  to  write  an  explicit 
expression  for  the  matrix  M.  It  is  easy  to  see  that  in  the  supersonic  case 

(49)  =  sign(«)  •  I 

where  /  is  a  3  x  3  unity  matrix. 

Some  algebra  reveals  that  in  the  subsonic  case  matrix  sign(A)  has  a  very  simple 
structure  as  well 


M^vh  _ 


sign(w)  0  0 

0  0  1/c 

0  c  0 


The  second  and  third  rows  mean  that  the  artificial  diffusion  added  to  the  momentum 
equations  is  proportional  to  the  fluctuation  of  the  pressure  equation  and  vice  versa. 

Remark  3.2.  The  auxiliary  variables  formulation  (37)  of  the  Euler  equations 
is,  perhaps,  the  simplest  way  to  write  the  Euler  equations.  The  usefulness  of  this 
formulation  will  become  even  more  evident  in  §§.^,6. 

4.  Construction  of  the  two-dimensional  Euler  scheme.  The  Euler  equa¬ 
tions  of  gas  dynamics  in  two  dimensions  can  be  written 


where 


Ut  +  F(u)a:  -f-  G{u)y  =  o, 


F(„)=  ;  G(„)  = 

puv  ^  ' 

\  puH 


pv^  +  p 

pvH 
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where  the  enthalpy  H  is  defined  by 


e  +  p 

the  speed  of  sound 

(54) 

and  the  pressure 

(55)  P  =  (7-  I)(e-P^"— ■  ) 

The  quasilinear  non-conservative  formulation  of  the  Euler  system  in  auxiliary  variables 
can  be  introduced  in  two  dimensions  as  well 


St  +  USj:  +  VSy  =  0 

put  +  puu^r:  +  pVUy  +P^  -0 

^  ^  pVt  +  pUVj:  +  pVVy  -\-py  =  0 

Pt  +  'tiPx  +  vpy  +  pc^{Uj,  +  =  0 

where  ds  =  dp  —  ^ . 

Remark  4.1.  Note  that  the  entropy  (s)  evolution  is  subject  to  the  two-dimensional 
advection  equation^  which  is  locally  decoupled  from  the  rest  of  the  system. 

The  fluctuation  of  the  system  (51)  defined  over  the  triangle  T 

(57)  R=  J  J  ut  =  -  J  J {Fx  +  Gy)  dx  dy  =  -St  +  Gy 

where  F^^Gy  are  some  averaged  values  of  the  flux  derivatives  over  the  triangle  T. 

Our  construction  of  the  truly  two-dimensional  Euler  scheme  utilizes  the  two  di¬ 
mensional  extension  (Roe,  Struijs  and  Deconinck  [16])  of  the  Roe  conservative  lin¬ 
earization  for  one-dimensional  case.  Therefore,  following  the  procedure  developed  in 
[16],  we  assume  that  the  quantity  which  varies  linearly  over  an  element  is  the  ‘^param¬ 
eter  vector” 


(58)  m  =  y/p{l,u,  v,HY 

and  its  averaged  value  on  the  triangle  T  (as  illustrated  on  Fig. 2)  is  given  by  the 
foUowing 


(59) 

mo  +  m3  +  m4 

Roe-averaged  quantities 

can  be  introduced 

(60) 

u  —  m2fmi 

V  =  m3 /mi 

H  =  Th^lfhi 

and 

(61) 

=  (7  -  1)[^  -  +  «^)] 
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Fluctuations  of  the  Euler  system  in  the  auxiliary  variables  can  be  presented  as 
(62)  r 


where 


ry  =  -StB  ■ 


u 

0 

0 

0  \ 

(  V 

0 

0 

0 

0 

u 

0 

1  5 

0 

V 

0 

0 

0 

0 

u 

0  ’ 

0 

0 

V 

1 

0 

C2 

0 

u  / 

0 

0 

C2 

V 

and  St  =  h?  12  is  the  area  of  the  triangle  T,  and 

(63)  px  —  2Tfl\(^7Tt\'^x 

(64)  ^  =  mi(m2)x  -  m2(Tni)x 

(65)  ^  =  rhi(m3)x  -  ni3{mi)x 

T  —  1 

(66)  px  =  - [(ih4(mi)x  +  mi{m4)x) 
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[{m4{mx)x  +  mi{m4)x)  +  (m2(”i2)a;  +  "i3(wi3)i:)] 


The  corresponding  terms  involving  derivatives  in  the  y  direction  can  be  written  in  the 
analogous  manner. 

Introducing  the  matrix 

/  1  0  0  l/c2  \ 

(67)  a  =  ®  i  “  % 

^  ’  V  0  1  vjc^ 

^^{u^  +  v‘^)|2  u  V  1/(7  —  1)  +  (it^  +  t;^)/(2c^)  j 


we  can  define 

(68) 

It  can  be  easily  verified  that 

(69) 


B?  = 

By  =  Cary 


By  = 

By  =  -S'^Gy, 


where  Fx,Gy  are  the  same  averaged  flux  derivative  values  as  defined  in  [16].  It  is  also 
obvious  that  the  entire  fluctuation 

(70)  B=B^  +  By  =  Ca(r^  +  ry)  = 

Consider  triangle  T  as  illustrated  in  Fig.2.  The  fluctuation  is  distributed  according 
to  the  following  formulae: 

=5uS  +^Ca(r^-f^) 

(71)  Sul+^  ^SuS  +^Cal(r-  +  r^)  +  (ry-ry)] 

5<+^  =5<  +^Ca(ry  +  ry) 
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Defining 


=  sign(i)r® 

^  =  sign(5)r^ 

we  obtain  the  scheme  that  is  similar  to  the  standard  Roe  dimensionally  split  scheme. 
The  only  difference  is  in  the  linearization  procedure. 

We  can  construct  now  a  (linearity  preserving)  second  order  accurate  scheme.  First, 
we  shall  introduce  vectors  with  their  elements  defined  by 

rf  =  rf  +  '^(qi)rf 
rf  =  rf  + 

for  i  =  1,2, 3, 4,  where 
(74) 

and  is  a  (non-compressive)  limiter. 

Substituting  for  r^,  in  (71)  and  (72)  we  obtain  a  genuinely  two-dimensional 

scheme,  which  is  also  linearity  preserving  (second  order  accurate  in  this  case).  This  is 
because  if 


r  =  0 


then 


as  well.  Therefore,  the  update  of  any  variable  in  all  the  three  vertices  of  the  triangle 
T  due  to  the  fluctuation  on  this  triangle  is  zero. 

The  resulting  scheme  is  also  conservative  since 

Some  attributes  and  properties  of  the  genuinely  multidimensional  schemes  wiD  be 
discussed  later  in  §9.3.  It  is  important  for  the  purpose  of  this  discussion  to  write  down 
the  explicit  expressions  for  the  matrices  sign(i4),  sign(jB).  Denote 


Mj:  =  sign(i) 
My  =  sign(5) 


For  matrix  the  distinction  should  be  made  between  two  cases 


(75) 

and  similarly 

(76) 


My  = 


Msub^ 

if  |it|  <  c 

if  |it|  >  c. 

if  |v|  <  c 

M^up^ 

if  |t;|  >  c, 
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where 


Mr'’  =  sigii({i)/, 


(78)  Mr'’  =  sign(t;)/. 

and  /  is  the  4x4  unity  matrix.  These  matrices  for  the  subsonic  case  appear  to  be 
very  simple  as  well 


Mr*”  = 


sign(tt)  0  0  0 

0  0  0  1/c 

0  0  sign(u)  0 

0  c  0  0 


Mr*”  = 


sign(u)  0  0  0 

0  sign(u)  0  0 

0  0  0  1/c 

0  0  c  0 


and  their  structure  may  have  a  very  important  meaning  (see  §9.3). 

5.  General  unstructured  grid  formulation.  Consider  an  unstructured  tri¬ 
angular  grid  covering  our  domain  and  triangle  T  belonging  to  this  grid  (see  Fig.4). 
Denote  by  hi  the  length  of  the  ith  face  of  this  triangle,  €{  a  unit  vector  along  the  ith 
face  in  the  clockwise  direction.  Denote  also  hi  a  unit  inward  normal  to  the  ith  face 
and  St  -  area  of  the  triangle  T.  In  this  section  we  shall  describe  first  the  general 
(unstructured  grid)  version  of  the  fluctuation-splitting  advection  scheme.  Then  we 
shall  construct  the  genuinely  two-dimensional  numerical  scheme  for  the  compressible 
Euler  system,  also  formulated  for  unstructured  meshes. 

5.1.  Advection  scheme.  Consider  a  linear  constant  coefficient  advection  equa- 


(81)  ut  +  \-  Vu  =  0 

Consider  a  non-orthogonal  coordinate  system  ^,7/  aligned  with  the  vectors  ei,e2  (see 
Fig.4)  and  introduce  the  following  quantities 

^0O^  A  •  7i2  ^  X  •  Ui 

(82)  “  =  — >  ^  = 


where 


6  =  ei  ■  h2  =  €2  ■  hi. 


Then  Eq.(81)  can  be  written  in  the  coordinate  system  ^,7/ 


Ut  -f  -h  —  0 


The  fluctuation  of  the  Eq.(81) 


R  =  R^-V 
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where 


=  -ST[ot{u3  -  ui)/hi] 

=  -5t[/3(u2  -  W3)//i2] 


The  fluctuation  distribution  formulae 

=  5i<  +  1{R^  -  R^) 

(87)  52^5+^  =  52k5  +  ^{R^  +  IR>) 

53U^+^  =  53^5  +  f  [(i?«  +  R^)  +  {R^  -  RP)] 

with  the  artificial  viscosity  terms  defined  as  follows 

.  R^  =  sign(a)i?^ 

^  R^  =  sigii(/3)i?’> 

result  in  a  positive  type  upwind  scheme.  Introduce  the  following  quantity 


Q  =  -- 


Define 


R^'  =  R^  +  R^^iQ) 

R^' =  R^  +  R^'^ 

and  substitute  R^’ ,  RP’  instead  of  R^,RP  in  (87), (88).  Using  the  identity 


RP^{Q)  =  -R^ 


we  can  rewrite  R^* ,  R^*  in  the  following  form 

(92) 

^  ^  R^'  ^  R^{\  -  ^'(g)). 

It  is  easy  to  verify  that  the  resulting  scheme  is  positive  if 

(93)  0<«’(g)<l,  0<^^<1 

Q 


and  linearity  preserving  if 


Note  also  that 


^(i)  =  i 


(95)  R  =  R^^  +  RP'  =R^  +  RP. 

Therefore,  the  constructed  scheme  is  conservative. 

There  are  two  differences  between  the  formulated  scheme  and  the  typical  fluctu¬ 
ation  splitting  advection  schemes  (see  [6], [22], [19]) 
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1.  no  distinction  is  made  between  triangles  with  two  outflow  faces  (Type  I)  and 
one  outflow  face  (Type  II) 

For  the  advection  problem  this  implies  more  computational  work  since  the 
limiter  function  has  to  be  evaluated  on  aU  the  triangles  (not  only  on  those  of 
Type  II).  However,  for  the  Euler  system  no  significant  savings  can  be  made 
by  this  distinction,  since  limiters  have  to  be  used  on  all  the  triangles  anyway 
(see  §5.2). 

2.  the  derivatives  are  approximated  by  finite  differences  along  the  faces  which 
were  chosen  arbitrarily  and  are  not  necessarily  both  inflow  or  outflow  (a  and 
/3  are  not  necessarily  of  the  same  sign). 

In  the  case  of  advection,  the  choice  of  faces  which  are  both  either  inflow  or 
outflow  for  derivative  approximation  provides  better  resolution  of  disconti- 
nidties.  For  systems,  however,  the  question  of  the  optimal  choice  is  more 
complex  (see  remark  in  §8).  Therefore,  we  presented  here  the  construction  of 
the  advection  scheme  with  no  assumption  made  about  the  signs  of  a  and 

5.2.  Euler  system.  The  auxiliary  variable  formulation  of  the  Euler  system  can 
be  written  in  the  following  vector  form 

+  17  •  Vs  =  0 

(96)  pUt  +J>U  •  VC7  +  Vp  =  0 

Pt  +  U  ■'Vp  +  pc^V  -U  =  Q 

where  U  —  (u,  v) 

Consider  agcdn  triangle  T  on  Fig.4  and  the  non-orthogonal  coordinate  system 
(^,p)  aligning  with  two  of  the  faces  of  this  triangle  (or  the  unit  vectors  61,62)  Define 
the  following  quantities 


(97) 

U  =  U-e2,  V  =  U-ei. 

and 

(98) 

U  ■  ill  „  V  -n^ 

We  can  rewrite  the  Euler  system  (96)  in  this  new  non-orthogonal  coordinates  (^,  rj) 

(99) 

St  -H  CLS^ 

pUt  +  poLi^  -1-  p(3Ur]  -b  =  0 

pVt  +  paV^  +  pfiVr,  +  p,,  =  0 

Pt  -b  ap^  -b  /3p„  -b  pc^{oL(^  -b  ^„)  =  0 

In  order  to  compute  the  fluctuation  of  the  system 

(100) 

r  =  -b  T*’' 

on  a  general  triangle 
presented  in  [16]  and 
rameter  vector 

T  we  foUow  again  the  two-dimensional  linearization  procedure 
assume  the  quantity  varying  linearly  on  the  triangle  is  the  pa- 

(101) 

m  =  y/p{l,u,v,H)'^ 
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and  its  averaged  value  on  the  triangle  T  is  given  by  the  following 

mi  +  m2  +  77X3 


Then  (similarly  to  the  structured  grid  case) 


where 

(104) 

(105) 

(106) 


—St 


=  2mi(mi)^ 


7  -  1 


as^ 

apU^  +  Vl 
apV^ 

+  c^p^l 


[(m4(mi)^  +  +  (7^2(7722)^  +  7773(7773)^)] 


p^-c  p^ 


pU^^ 


_  f  Ph\^(  ^1(^2)^-7772(7711)^ 


7^1(7713)^  -  7n3(77Zi)^ 


can  be  evaluated  in  a  similar  way.  Defining  the  following  set  of  conservative  vari¬ 
ables 


and  denoting  R  the  fluctuations  of  the  Euler  system  in  these  variables  on  triangle  T, 
it  can  be  easily  verified  that  the  fluctuations  in  the  following  relation  holds 


(109) 

where 


(110) 


R  =  Ca-r 


1  0  0  l/c2 

ZY  1  0  ZY/c2 

V  0  1  V/c2 

f72/2  a  ^  1/(7- l)  +  t/V(2c2) 


The  fluctuation  distribution  formulae 


(111) 


5iu”+'  =Siu^  +iC„(r«-f«) 

S2U^-^^  =S2U^  +^Ca{r^  +  f^) 

=53*5  +ia[(r«+rO  +  (^"-^")] 


will  result  in  an  upwind  scheme  provided  we  define  the  artificial  viscosity  terms 
However,  before  doing  this  we  would  like  to  rewrite  (111)  for  the  update  of  the  ‘‘reg¬ 
ular”  conservative  variables  u  =  (/>, />u,/97;,e)^.  Note  that 


u  =  0  • 


(112) 
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where 


10  0  0 
0  niie  ni/e  0 

0  nf  /  ^  ^2  /  ^  0 
0  0  0  1 


Introducing  the  following  matrix 


'Ci  =  Q-Ca  = 


1 

0 

0 

l/c^ 

u 

nf/e 

n\/e 

ujc^ 

V 

n^/e 

nl/0 

n/c^ 

U^/2 

a 

5 

1/(7-1)  +  C7V(2c2) 

we  obtain  the  fluctuation  distribution  formulae  for  the  “regular”  conservative  variables 

(115)  52U^'  + 

53U^+'  =531*?  +iCf[(r«  +  r«)  +  (r’'-r^)] 


The  relationship  between  the  artificial  viscosity  and  the  fluctuations  is  given  by 


(116) 

where 

(117) 
and 

(118) 


=  Mr,T^ 


if  |d|  <  c 
if  |d|  >  c, 


Mr,  = 


if  |/3|  <  C 

M-P,  if|5|>c, 


Msub  ^ 


sign(d)  0  0  0 

0  0  0  1/c 

0  0  sign(d)  0 

0  c  0  0 


and 

(121) 


Msub  _ 
V 


sign(;3)  0  0  0  ^ 

0  sign(5)  0  0 

0  0  0  1/c  ’ 

0  0  c  0  ; 


=  sign(a)/, 
=  sign(5)/. 


A  linearity  preserving  genuinely  two-dimensional  scheme  can  be  constructed  by 

with  their  elements  defined  by 


introducing  the  vectors 

(123) 

for  i  =  1 5  2, 3, 4,  where 

(124) 


=  r>  +  ^'(g,)^!' 


V  V  I 

r-  =  r-  + 


9i 


(li-  — V 
rj 


and  substituting  them  instead  of  respectively  in  (115), (116). 

Remark  5.1.  In  order  to  construct  a  non- oscillatory  linearity  preserving  scheme 
for  the  Euler  equations^  derivatives  can  be  approximated  along  any  two  out  of  three 
faces  of  the  triangle.  However ^  the  resolution  of  such  features  of  the  flow  as  shocks, 
contact  discontinuities  and  shear  layers  can  be  affected  by  the  particular  choice  (see 
^8).  Shocks  are  being  resolved  better  (i.e.  are  represented  by  sharper  numerical  layers) 
if  the  numerical  derivatives  are  computed  along  those  two  faces  of  the  triangle  which 
are  the  closest  to  the  direction  of  the  shock.  The  same  is  true  for  the  shear  layers  and 
contact  discontinuities,  although  the  two  faces  closest  to  the  discontinuity  direction  in 
this  case  will  also  be  either  two  inflow  or  outflow  faces  of  the  triangle. 


6.  Outline  of  the  extension  to  three  dimensions.  A  detailed  description  of 
the  three-dimensional  schemes  both  for  the  scalar  advection  and  Euler  system  will  be 
given  elsewhere.  In  this  paper  we  only  outline  briefly  their  construction  in  the  case  of 
structured  Cartesian  grids.  Each  cube  having  grid  nodes  as  its  vertices  can  be  divided 
into  two  prisms.  The  prisms  can  be  divided  into  three  tetrahedra  each.  Any  of  the 
obtained  tetrahedra  will  have  3  of  its  edges  parallel  to  the  x,  y  and  J2r  axes  respectively. 
We  shall  utilize  this  in  this  section  to  make  the  presentation  free  of  unnecessary  details 
related  to  the  non-orthogonal  coordinate  systems.  The  construction  of  the  fluctuation¬ 
splitting  schemes  both  for  the  scalar  advection  and  the  Euler  system  will  be  outlined 
for  the  case  of  the  tetrahedron  T  as  illustrated  on  Fig. 5. 

6.1.  Advection  scheme  in  three  dimensions.  Consider  a  scalar  constant  co¬ 
efficient  advection  equation  in  three  dimensions 


(125) 

Ut  +  +  bUy  +  cuz  =  0 

The  fluctuation  of  this  equation  on  the  tetrahedron  T  (see  Fig. 5),  whose  volume  is 
V  =  h^/6. 

(126) 

R  =  +  Ry  +  R^ 

where 

(127) 

R^  =  -^[a{u3  -  U4)] 

-  ''is)] 

R^  =  -  U2)] 

The  fluctuation  distribution  formulae  are  given  as  follows 

=  Vu^  +zc,{R^  +  R^) 

=  Vu^  +LC^[{Ry  +  Ry)+{R^^R-)] 

=  +rc^[{R^  +  R^)  +  {Ry -Ry)] 

=^Vu2  +^Ca{R^-R^) 
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(128) 


(130c) 

where 

(131) 
and 

(132) 
or 


mj  +  m: 

Rz  +  +  [Ry]t 


f  jR%  if  sign(i?‘)  =  sign(f2'*), 

1  0,  otherwise, 


0,  if  sign(B‘)  =  sign(jR''), 
—R\  otherwise. 


(133)  =  [R^]+  -  W 

where  i  and  k  stand  for  one  of  (though  different)  x^y  or  z  each. 

Then  we  define  the  following  quantities 

=R^  +  ^{Qy)Ry  +  ^{Q^)R^ 

(134)  Ry"  ^Ry  +  ^Iq^)R^  +  ^Iq^)R^ 

RZ*  =  Rz  +  $(Q^)i2^  +  ^lQy)Ry 

where  ^  is  a  non-compressive  limiter. 

Substituting  R^*^Ry* ,R^*  into  (128), (129)  instead  of  R^^Ry^R^  respectively  we 
obtain  a  linearity  preserving  (second  order  accurate  in  this  case)  positive  advection 
scheme. 

To  demonstrate  this  we  can  assume  without  loss  of  generality  that 


(135) 

Qy-Q^>Q 

Q^{Qy  +  Q^)<0 

Then 

(136) 

Qy  =  -R^/(Ry  +  R^) 
=  ilQy 

Q^  =  Qy 
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’^{Qy)Ry  +  ^iQ^)R^  =  -^{^)R^ 

^{Q^)R^  =  -^{Qy){Ry  +  R^) 

R-'  =  (1  -  ^{^))R^ 

Ry*  =  (1  -  if!{Qy))Ry 
R^'  =  (1  - 

It  can  be  concluded  that  if  is  a  non- compressive  limiter,  then  the  constructed  scheme 
is  positive. 

If  the  fluctuation  R  vanishes  then 

(139)  =  R^  +  R^ 
and,  therefore, 

(140)  Q-^Qy=.Q-  =  L 
Then  it  is  easy  to  see  that 

(141)  R^*  =  R^*  =  R^*  =  0, 

provided  ^^(1)  =  1,  i.e.  the  constructed  scheme  is  linearity  preserving. 

6. 2,  Euler  system.  Euler  equations  in  three  dimensions 

(142)  Ut  -h  F{u)j:  -f  G{u)y  +  if(u)^  =  o 
where 


In  this  case 

(137) 

and  therefore 

(138) 


F(u)  = 


(  P  \ 
pu 


u  = 

pv 

pw 

\  ^  ) 

(  pu  \ 

(  pv 

(  pw  \ 

pv?  +  p 

puv 

pwu 

puv 

;  G(«)  = 

pv^  -|-  p 

;  H{u)  = 

pwv 

puw 

pwv 

pw^  +  p 

u(e  +  p)  ) 
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Introducing  the  auxiliary  variables  (5,  t/,  u,  t/;,p)  we  can  rewrite  these  equations  in  the 
following  non- conservative  form 

St  +  -h  VSy  d-  WSz  =  0 

pUt  +  pUUa:  +  pVUy  -h  pWU;r  +  =  0 

(143)  pVt  +  pUV:c  +  pVVy  +  pWVz  +Py  =  0 

pWt  +  pUWa;  +  pVWy  +  pWWz  +p^  =  0 

Pt  +  Upa:  +  VPy  +  WPz  +  pC^{U:c  +  =  0 
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Fltictiiatioii  on  tetrahedron  T  (see  Fig.5)  is 


(144) 


-III  Ut  =  ~  j  J  J +  Hz)dx  dy  dz 

=  -VriK  +  ^y  +  Hl] 


Flnctnation  in  auxiliary  variables  on  each  tetrahedron  can  be  represented  as  a  sum  of 
three  parts 


(145) 


r  =  r®  +  +  T*-^ 
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1/(7  -  1)  +  UV(2c) 

The  matrix  Ca  relating  the  fluctuations  in  conservative  variables  R  to  those  in  the 
auxiliary  variables 


(146) 


It  was  pointed  out  in  [16]  that  the  two-dimensional  conservative  linearization  presented 
there  extends  to  three  dimensions  in  a  straightforward  way.  Again,  we  assume  that 
the  quantity  that  varies  linearly  on  the  tetrahedron  T  is  the  “parameter-vector”  and 
the  averaged  quantities  can  be  defined  in  a  manner  similar  to  the  two-dimensional 
case.  The  fluctuation  in  the  auxiliary  variables  on  T  can  be  written  as  follows 

(147)  r  =  -St  •  a  •  r""  =  -StA  • 

The  fluctuation  distribution  formulae  in  this  case  are 


(148) 


=Vu^  +fCa(r^  +  r") 

=Vu^  +^Ca[iry  +  fy)  +  {r^  -  r^)] 
=Vu'^  +^Ca[ir^  +  r^) +  iry  -  fy)] 
=Vu2  +^Ca{r^-f^) 


Similarly  to  the  two-dimensional  case,  the  relationship  between  the  the  artificial  vis¬ 
cosity  needed  in  x-direction  to  obtain  the  upwind  scheme  and  is  the  following 

(149)  =  sign(A)r^. 

and  can  be  evaluated  in  a  similar  way. 

Denoting 


it  can  be  verified  that 
(150) 


M^  = 


Mx  =  sign(i) 


if  |u|  <  c 
M-P,  if  |u|  >  c, 
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(155) 


[^f]j  +  [r\]: 

r^i  +  [rf]t  +  [r\]t 


where  [. .  .]J,  [. .  [. .  .]*  are  defined  by  (131), (132).  Now  we  can  introduce 


rf  =  7-f  +  «’(gf)rf  +  ^{ql)rl 

(156)  rf  =  rf  +  $(gf)rf  +  ^(gf  )rf 

rf  =  rf  +  5’(gf)rf  +  ^(gf>f 

for  i  =  1, 

Substituting  ,r^  ,r^  instead  of  into  the  fluctuation  distribution  for¬ 

mulae  (148)  and  (149)  we  obtain  a  linearity  preserving  (second  order  accurate  for  the 
case  of  structured  meshes)  genuinely  three-dimensional  scheme  for  the  Euler  equations. 
Note  that 


(157) 


r  =  +  r^*  +  . 


Therefore,  the  constructed  scheme  is  conservative.  It  is  also  easy  to  see  that  it  is 
linearity  preserving. 

7.  Multigrid  solver.  One  of  the  most  attractive  properties  of  the  constructed 
multidimensional  scheme  for  the  Euler  equations  is  that  the  Collective  Gauss-Seidel 
relaxation  is  stable  when  applied  directly  to  the  high-resolution  discrete  equations. 
This  can  be  utilized  to  construct  a  very  simple  and  efficient  steady-state  multigrid 
solver.  In  this  work  we  would  like  just  to  illustrate  the  basic  advantages  of  the  multigrid 
solver  that  uses  the  constructed  multidimensional  scheme.  Therefore,  we  do  not  discuss 
here  any  issues  concerning  the  multigrid  algorithm  for  general  unstructured  meshes, 
but  address  the  interested  reader  to  the  literature  (see,  for  instance,  [11]). 

We  shall  present  briefly  in  this  section  the  simple  multigrid  algorithm  used  in  the 
numerical  experiments  presented  in  this  work. 
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Description  of  the  algorithm.  The  multigrid  cycle  used  in  this  work  is  W-cycle 
and  the  entire  algorithm  can  be  implemented  either  as  Cycling  or  as  Full  Multigrid 

{FMOy 

Relaxation,  The  algorithm  employs  the  Collective  Gauss-Seidel  relaxation  as  a 
smoother.  In  order  to  update  the  solution  at  each  node  a  system  of  four  nonlinear 
equations  has  to  be  solved.  This  can  be  done  using  the  Newton  method.  One  iteration 
at  each  node  is  sufficient  for  this  purpose  (see  [18]  for  the  discussion  on  this  matter). 
The  ordering  of  the  relaxation  can  be  Red-Black,  lexicographic  etc. 

Restriction  and  prolongation.  Assume  that  we  have  a  hierarchical  sequence  of  tri¬ 
angular  grids,  i.e.  that  any  triangle  on  the  coarser  grid  is  a  union  of  four  triangles 
on  the  finer  grid.  The  natural  choice  for  prolongation  in  this  case  would  be  a  linear 
interpolation  along  each  face  of  the  triangle  belonging  to  the  coarser  grid.  The  re¬ 
striction  in  this  case  can  be  just  combining  the  fluctuations  from  the  four  finer  grid 
triangles  forming  the  coarse  grid  one.  However,  we  consider  only  Cartesian  grids  in 
the  numerical  experiments  presented  in  this  work  (see  Figs.2,3).  In  this  case  we  can 
use  more  conventional  prolongation  and  restriction  operators:  bilinear  correction  in¬ 
terpolation  and  FuU- Weighting  of  the  residuals  (see  [1]).  Term  “residual”  is  used  here 
in  its  usual  sense:  residual  of  the  discrete  equation  at  each  gridpoint  (constructed  from 
the  contributions  from  the  triangles  having  this  gridpoint  as  a  common  vertex). 

Previous  relevant  results.  Application  of  multigrid  methods  for  the  advection  equa¬ 
tion  was  studied  in  detail  in  [17], [18]  in  conjunction  with  the  genuinely  two-dimensional 
control- volume  type  scheme  developed  there.  It  was  demonstrated  there  that  applying 
the  2FMG  —  kF(2, 1)  algorithm  is  sufficient  to  obtain  high  quality  solutions  (i.e.  solu¬ 
tions  with  sharply  resolved  discontinuities  and  second  order  accuracy  both  in  smooth 
regions  and  in  discontinuity  location)  for  the  advection  equation.  Many  of  the  conclu¬ 
sions  of  [17], [18]  apply  directly  to  the  Euler  system  case  considered  in  this  work. 

Efficiency  of  the  multigrid  solver  for  advection-dominated  problems.  It  was  shown 
in  [1]  that  the  residual  reduction  per  multigrid  cycle  for  an  advection  (dominated) 
problem  cannot  in  general  be  better  than  .75  for  the  second  order  accurate  discretiza¬ 
tion  (.5  for  the  first  order).  This  is  because  the  coarse  grid  provides  only  a  fraction 
(.75  or  .5  respectively)  of  the  needed  correction  for  some  components.  The  ability  of 
the  multigrid  algorithm  to  achieve  a  rate  of  convergence  close  to  .75  (for  the  second 
order  accurate  approximation)  means  that  the  obstacle  towards  achieving  a  better 
efficiency  are  not  the  smoothing  properties  of  the  relaxation  but  rather  an  insufficient 
coarse  grid  correction  for  certain  components,  (see  §9.2.2  for  further  discussion  on  this 
issue). 

8.  Numerical  experiments.  The  purpose  of  the  numerical  experiments  re¬ 
ported  in  this  section  is  to  verify  the  robustness  of  the  constructed  scheme  and  the 
quality  of  the  numerical  solutions  obtained  by  its  means.  Some  preliminary  experi¬ 
ments  illustrating  the  performance  of  the  multigrid  algorithm  uzing  this  scheme  are 
presented  as  well. 

8.1.  Supersonic  flow  in  a  channel  with  a  bump.  The  test  case  considered 
here  is  a  supersonic  (Mach=2.9)  flow  in  channel  with  a  circular  bump.  The  bump  is 
located  at  the  lower  wall  of  the  channel  at  1  <  x  <  2.  and  its  surface  is  a  circular  arch 
of  7r/3  and  radius  1.  Note  that  the  actual  shape  of  the  domain  is  a  rectangle.  The 
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influence  of  the  bump  on  the  flow  is  imposed  through  the  boundary  conditions:  the 
velocity  component  normal  to  the  surface  of  the  bump  at  a  certain  location  is  being 
reflected. 

The  first  experiment  uses  a  grid  of  size  200  x  40  points.  The  density  contour  plots 
of  the  steady-state  solution  are  presented  on  Fig. 6(a).  The  scheme  used  is  the  one 
given  by  (71), (72),  (73)  in  §4  with  the  minmod  limiter. 

The  second  experiment  presented  in  Fig. 6(b)  corresponds  to  the  same  settings, 
except  that  the  grid  is  twice  finer  (400  x  80  points).  As  it  is  expected,  the  grid 
refinement  results  in  a  better  resolution  of  the  flow  features. 

The  third  experiment  (Fig.6(c))  is  performed  on  the  grid  of  the  same  size  as 
the  first  one.  However,  the  triangulation  is  as  illustrated  on  Fig. 3  and  the  scheme 
employed  is  constructed  according  to  (115),(116),(123)  in  §5  and  the  derivatives  are 
approximated  using  the  diagonal  and  parallel  to  the  x  direction  faces  of  the  triangles. 
Even  though  the  grid  is  structured,  the  most  essenticd  feature  of  the  unstructured  grid 
formulation  of  the  scheme  -  use  of  the  non-orthogonal  coordinate  system  is  present. 

The  purpose  of  this  experiment  is  to  test  the  performance  of  such  a  scheme  and  to 
illustrate  the  effect  of  the  scheme  and  triangulation  on  the  resolution  of  discontinuities. 
It  can  be  seen  in  Fig. 6(c)  that  the  shock  reflected  from  the  upper  wall  is  resolved  better 
than  the  stronger  one  that  is  incident  to  the  upper  wall.  This  is  because  for  this 
scheme/grid  combination,  discontinuities  whose  direction  is  close  to  the  grid  diagonal 
(upper-left  to  lower-right)  will  be  resolved  very  weU  (in  addition  to  those  which  are 
close  to  the  horizontal  or  vertical  directions). 

8.2.  Transonic  flow  over  a  circular  bump.  The  testcase  considered  here  is 
a  transonic  flow  (free-stream  Mach=  .9)  over  a  flat  wall  with  a  bump  (Fig. 7).  The 
surface  of  the  bump  is  a  circular  arch  of  7r/3  and  radius  1  and  its  location  is  between 
3.5  <  X  <  4.5.  Again,  in  order  to  keep  the  experiments  simple  at  this  stage  of  work, 
the  bump  is  treated  the  same  way  as  in  the  previous  experiments.  The  grid  is  200  x  200 
points  and  the  scheme  used  is  presented  in  §4.  The  shock  of  the  “fish-tail”  shape  can 
be  clearly  observed  in  Fig. 7. 

8.3.  Low  Mach  number  flow  over  a  circular  bump.  Here  we  present  a 
numerical  experiment  concerning  a  low  Mach  number  (=.l)  flow  over  a  flat  waU  with 
a  circular  (arch  of  7r/3  and  radius  2)  bump.  Here  as  well  as  in  the  previous  case 
the  presence  of  the  bump  is  imitated  through  the  appropriate  boundary  conditions. 
The  grid  is  200  x  200  points.  The  density  contours  of  the  steady-state  solution  are 
presented  in  Fig. 8. 

8.4.  Multigrid  algorithm.  To  illustrate  the  performance  of  the  multigrid  algo¬ 
rithm  we  consider  here  the  well  known  testcase  of  a  shock  reflecting  from  a  flat  wall. 
The  multigrid  algorithm  involves  five  grids  (levels):  the  finest  consists  of  129  x  33 
points,  the  coarsest  is  9  x  3  points. 

The  multigrid  algorithm  is  based  on  the  scheme  presented  in  §4  used  with  the 
lexicographic  Gauss-Seidel  relaxation.  The  restriction  and  prolongation  procedures 
are  the  standard  Full  Weighting  of  the  residuals  and  bilinear  correction  interpolation. 
The  first  experiment  performs  one  W(2, 1)  cycle.  Fig.9(a)  presents  the  initial  guess 
and  Fig. 9(b)  -  the  numerical  solution  obtained  by  one  W{2,1)  cycle.  It  is  remarkable 
that  such  little  computational  work  is  sufficient  to  make  the  reflected  shock  clearly 
visible.  The  numerical  solution  to  this  problem  obtained  by  the  2FMG  -  W{2,1) 
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algorithm  is  presented  on  Fig.9(c). 

Note  that  in  this  case  the  flow  is  aligned  with  the  x-direction  in  a  significant  part 
of  the  domain.  In  this  case  the  artificial  viscosity  in  the  cross-stream  direction  in  the 
entropy  and  ti-momentnm  equations  (see  §4)  vanishes.  Therefore,  no  smoothing  can 
be  obtained  in  the  ^-direction  in  some  components.  A  multigrid  algorithm  utilizing 
the  time-stepping  type  relaxation  can  deal  with  such  a  situation  only  using  the  semi¬ 
coarsening  technique.  Our  algorithm  employs  the  Gauss-Seidel  relaxation.  Therefore, 
it  offers  a  much  simpler  and  more  efficient  treatment  of  this  problem:  relaxation  with 
lexicographic  ordering  in  the  stream  direction. 

The  rate  of  convergence  observed  in  this  testcase  as  weU  as  in  other  simple  ex¬ 
periments  concerning  a  variety  of  flow  regimes  is  very  close  to  .75  (see  §§7, 9.2.2  for  a 
discussion). 

9.  Discussion  and  conclusions. 

9.1.  Summary  of  the  current  work.  A  new  two-dimensional  high-resolution 
(at  the  steady-state)  scheme  for  the  compressible  Euler  equations  was  presented.  It 
is  triangle-based  and  can  be  formulated  with  the  same  degree  of  simplicity  both  on 
structured  and  unstructured  grids.  The  main  advantage  of  this  scheme  is  that  Gauss- 
Seidel  relaxation  can  be  applied  directly  to  the  resulting  discrete  equations.  This 
allows  construction  of  a  simple  and  efficient  multigrid  steady-state  solver. 

A  remarkable  property  of  the  constructed  scheme  is  also  its  very  compact  stencil: 
it  involves  only  the  immediate  neighbors  of  the  point  of  interest. 

A  variety  of  flow  regimes  (supersonic,  transonic  and  low  Mach  number  flow)  were 
considered  in  the  numerical  experiments  to  verify  the  quality  of  the  solutions  obtained 
by  means  of  the  new  scheme  and  to  demonstrate  the  efficiency  of  the  multigrid  algo¬ 
rithm. 

Generalization  of  this  scheme  to  three  dimensional  tetrahedral  meshes  is  presented 
briefly  as  well. 

9.2.  Future  work. 

9.2.1.  Compressible  Navier-Stokes  equations.  The  extension  of  this  scheme 
to  the  case  of  the  compressible  Navier-Stokes  equations  is  straightforward.  It  is  inter¬ 
esting  to  note  that  the  artificial  viscosity  present  in  the  momentum  equations  in  the 
subsonic  case  contains  already  terms  of  the  type  {ux  +  Vy)x  and  {ux  +  Vy)y  that  appear 
in  the  physical  viscosity  of  the  compressible  Navier-Stokes  equations. 

9.2.2.  Further  improvement  of  the  multigrid  efficiency.  As  it  was  men¬ 
tioned  in  §7  the  main  obstacle  towards  the  further  improvement  of  the  multigrid 
efficiency  is  the  following  fact:  for  the  hyperbolic  problems  the  coarse  grid  correction 
is  not  sufficient  for  certain  error  components. 

This  difficulty  was  already  addressed  in  the  literature  and  some  techniques  to 
improve  the  multigrid  efficiency  were  developed  in  [2].  Therefore,  one  possibility  is  to 
adapt  these  techniques  for  our  case  -  compressible  Euler  equations. 

We  shall  mention  here  another  way  to  deal  with  this  difficulty.  The  steady-state 
compressible  Euler  equations  in  the  subsonic  case  are  very  sinoiilar  from  the  mathe¬ 
matical  point  of  view  to  the  incompressible  Euler  equations  -  both  problems  are  of 
the  mixed  eUiptic-hyperbolic  type. 
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There  exist  very  efficient  steady-state  multigrid  solvers  for  the  incompressible  flow 
equations,  which  rely  on  the  primitive  variables  formulation  (velocities  and  pressure). 
Most  of  them  employ  staggered  grid  discretization  of  the  equations.  A  typical  efficiency 
of  such  algorithms  for  a  low  Reynolds  number  case  is  comparable  to  that  for  the  Poisson 
equation.  However,  the  convergence  rate  deteriorates  for  the  advection  dominated 
(high  Reynolds  number)  flow.  Typical  numerical  schemes  for  incompressible  flow  are 
constructed  in  such  a  way  that  a  certain  “separation”  of  the  elliptic  and  hyperbolic 
(advection)  factors  of  the  equations  is  achieved  in  this  case  (see  [1]).  Assume  that  the 
flow  structure  is  simple  enough  (i.e.  there  is  no  recirculation,  etc.),  and  it  is  possible 
to  perform  relaxation  in  the  stream  direction.  Then  one  can  still  recover  the  optimal 
efficiency  efficiency  of  the  multigrid  solver  (see  [25]). 

However,  constructing  a  discretization  that  achieves  a  similar  separation  of  the 
advection  and  “full-potential”  factors  of  the  compressible  Eider  equations  appeared 
to  be  a  very  difficult  problem  despite  the  similarity  of  the  steady-state  equations  for 
compressible  and  incompressible  cases.  A  canonical  variable  formulation  of  the  Euler 
equations  for  both  cases  was  suggested  recently  by  Ta’asan  [23].  This  formulation 
facilitated  the  construction  of  the  discrete  scheme  (using  staggered  grids)  and  of  the 
relaxation  procedure  that  separate  a  treatment  of  the  advection  and  full-potential 
factors  (see  [24]).  The  resulting  multigrid  algorithm  achieves  a  very  good  efficiency 
(the  same  as  for  the  full-potential  equation)  in  the  subsonic  case  provided  it  is  possible 
to  perform  relaxation  in  the  flow  direction. 

One  of  the  future  possible  directions  is  to  modify  the  scheme  presented  in  this  work 
so  that  it  will  achieve  a  similar  efficiency  for  the  case  of  triangular  (non- staggered  and, 
possibly,  unstructured)  grids.  This  is  possible  because  some  very  important  similarities 
between  the  present  scheme  and  the  existing  discretizations  for  the  incompressible 
flow  equations  used  in  very  efficient  multigrid  solvers  (see,  for  instance,  [1])  can  be 
observed  (see  §9.3).  The  resulting  discretization  will  rely  on  the  primitive  rather  then 
the  canonical  variables.  The  latter  is  crucial  for  the  further  generalization  of  the  solver 
to  the  compressible  Navier-Stokes  equations. 

9.2.3.  Time-dependent  problems.  The  constructed  scheme  is  capable  of  pro¬ 
ducing  high-resolution  (second  order  accurate  on  the  structured  grids)  steady-state 
solutions  for  the  compressible  Euler  equations.  A  very  natural  way  to  extend  this 
capability  to  the  transient  problems  is  to  use  a  Lax-WendrofF  type  modification  of  the 
constructed  scheme,  i.e.  to  scale  the  artificial  viscosity  in  such  a  way  that  it  will  cancel 
the  time  error  of  the  forward  Euler  time-stepping  due  to  the  equation  (or  system  of 
equations)  being  solved. 

To  illustrate  this  on  the  simple  case  of  the  two-dimensional  advection  consider  the 
scheme  for  the  structured  grids  presented  in  §2.1.2.  The  use  of  the  following  artificial 
viscosity 


(158) 


=  laR^ 

Ry  =  i-hRy. 


instead  of  the  one  given  by  (5)  results  in  a  high-resolution  scheme  (constructed  by 
substituting  ,jR^  defined  by  (6)  instead  of  that  is  second  order  accurate 

in  time  and  space.  However,  it  is  easy  to  see  that  in  general  this  scheme  is  not  of  the 
positive  type  anymore  (see  also  [6]). 
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For  the  case  of  the  Euler  system,  modifying  the  construction  of  the  high-resolution 
scheme  presented  in  §4  by  substituting  the  artificial  viscosity  defined  as  follows 

=  IAt^ 

instead  of  (72),  results  in  a  Lax-Wendroff  type  scheme.  However,  the  applicability  of 
this  scheme  is  restricted  to  problems  with  no  strong  shocks. 

It  is  interesting  to  note  that  the  two-dimensional  schemes  by  Colella  [3], [4],  LeV- 
eque  [10],  Radvogin  [12]  can  be  interpreted  as  schemes  of  the  Lax-Wendroff  type. 
However,  the  nonlinear  mechanism  that  enables  them  to  deal  with  (strong)  shocks 
relies  on  one-dimensional  limiters  and  characteristic  variables  in  the  x  and  y  direction, 
which  introduces  an  essential  flavor  of  dimensional  splitting. 

The  first  step  towards  the  construction  of  a  genuinely  multidimensional  time- 
space  accurate  and  robust  Euler  scheme  should  be  to  construct  a  positive  time-space 
accurate  scalar  advection  scheme  of  the  type  Lax-Wendroff  type.  This  is  one  of  the 
directions  currently  being  investigated. 

9.3*  What  is  a  truly  multidimensional  scheme  ?.  Perhaps,  it  wiU  be  easier 
to  answer  this  question  after  the  nonosciUatory  time-space  accurate  multidimensional 
scheme  for  the  compressible  Euler  equations  is  constructed.  However,  some  attributes 
of  the  truly  multidimensional  schemes  can  be  seen  already  in  the  scheme  constructed 
here. 

One  of  the  main  characteristics  of  a  truly  multidimensional  scheme  for  fluid  dy¬ 
namics  is  that  a  multidimensional  operator  such  as  the  divergence  operator  (present 
in  the  pressure  equation,  the  auxiliary  variables  formulation  of  the  Euler  system  (56)) 
should  be  represented  in  a  difference  scheme  as  whole.  This  is  not  the  case  for  the 
dimensionally  split  schemes.  It  can  be  clearly  seen  in  the  example  of  the  dimensional 
upwinding  scheme  (71), (72)  that  terms  and  Vy  are  separated  (when  contributing 
to  the  artificial  viscosity  of  the  momentum  and  pressure  equations). 

Consider  the  case  of  subsonic  flow.  In  this  case  the  fluctuation  of  the  pressure 
equation  contributes  to  the  artificial  diffusion  of  the  momentum  equations.  Let  us 
modify  slightly  the  scheme  constructed  in  §4  and  compute  the  fluctuation  of  the  pres¬ 
sure  equation  contributing  to  the  artificial  viscosity  terms  of  the  momentum  equations 
according  to  the  following 

(160)  rf  =  r|  =  r J  -h  r|. 

The  only  difference  between  (160)  and  (73)  is  that  the  limiter  function  is  set  to  unity. 
Although 

rf  +  rf  =  2r4  ^  u, 

such  a  scheme  is  conservative,  since  (160)  is  used  only  in  the  artificial  viscosity.  This 
is  because  the  role  of  the  artificial  viscosity  can  be  interpreted  as  “redistribution”  of 
the  fluctuations  split  by  the  central  scheme  -  to  subtract  a  certain  amount  from  one 
vertex  of  the  triangle  and  to  add  the  same  amount  to  another  vertex.  This  does  not 
violate  conservation,  regardless  of  what  amount  has  been  “redistributed”. 

Recall,  that  for  some  very  popular  schemes  for  incompressible  flow  (see,  for  in¬ 
stance,  [1])  the  residual  of  the  continuity  equation  contributes  to  the  velocities  update 
in  such  a  way  that  the  discrete  vorticity  remains  unchanged. 
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The  “physical”  meaning  of  this  is  that  the  change  in  (evolution  of)  the  vorticity 
is  only  due  to  the  residuals  of  the  momentum  equations,  as  one  would  expect.  The 
“mathematical”  meaning  of  this  fact  is  that  the  elliptic  factor  of  the  incompressible 
Euler  equations  is  being  discretized  and  relaxed  as  such.  This  allowed  construction  of 
a  highly  efficient  multigrid  solver  (see  [1]). 

What  can  be  observed  clearly  in  the  subsonic  case  is  that  the  fluctuation  of  the 
pressure  equation  (being  a  part  of  artificial  diffusion  of  the  discretized  momentum 
equations)  contributes  to  the  change  of  the  velocities  (divergence  update)  according 
to  exactly  the  same  pattern  as  the  residual  of  the  continuity  equation  in  the  schemes 
for  discretizing  the  equations  of  incompressible  flow. 

Thus,  the  scheme  presented  in  this  work  establishes  a  link  between  some  standard 
techniques  used  for  solving  equations  of  the  incompressible  flow  and  the  class  of  upwind 
schemes  commonly  used  for  compressible  flow  computations. 
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Fig.  6.  Supersonic  flow  in  a  channel  over  a  circular  bump:  a)  grid  200  x  40  pis.,  triangulation  /, 
scheme  from  14;  h)  the  same,  except  the  grid  400  x  80  pts,;  c)  grid  200  x  40  pts.,  triangulation  II, 
scheme  according  to  §5,  which  relies  on  the  approximate  derivatives  along  the  nonorthogonal  faces  of 
the  triangles. 
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